clear;

%import data as exported by STARexportformatlab.do
regressors = csvread('regressors.csv',1,0); 
IFbounds= csvread('IFbounds.csv',1,0);  

n=length(regressors);

missingoriginal=regressors(:,5);
missingoriginal=logical(missingoriginal);
notmissingoriginal=~missingoriginal;


%%
%pointestimates
%generate "re-sampled" observations by taking the original sample
sample=1:n;
%run script to get point-estimates for bounds
calculatebounds;
boundsoriginal=bounds

%%
%replicate the analysis by drawing bootstrap samples
%store estimates for each replication in the matrix replicates
replications=1000;
replicates=zeros(replications,4)

for i=1:replications
    sample=ceil(rand(n,1)*n);
    calculatebounds;
    replicates(i,:)=bounds;
end

%variance of estimates across replications
variancebounds=cov(replicates);
stdbounds=std(replicates);
stdboundssigned= [stdbounds(1:2) -stdbounds(3:4)];

%%
%construct "uniform confidence set" as described in section 5.2 of the
%paper

%find delta and epsilon
alpha=0.05;
Zdraws=randn(5000,4) * mpower(variancebounds,0.5);

a=1;
delta=0;
while a>alpha
    delta=delta+0.0005;
    confdraws=Zdraws + delta* ones(5000,1)*stdboundssigned;
    er=confdraws(:,1) < 0 | confdraws(:,2)<0 | confdraws(:,3)>0 | confdraws(:,4)>0;
    a=mean(er);
end
%bounds for a superset, containing the identified set with probability
%1-alpha
G2hatbounds=boundsoriginal + delta*stdboundssigned


a=1;
epsilon=0;
while a>alpha
    epsilon=epsilon+0.0005;
    confdraws=Zdraws - epsilon* ones(5000,1)*stdboundssigned;
    er=confdraws(:,1) > 0 | confdraws(:,2)>0 | confdraws(:,3)<0 | confdraws(:,4)<0;
    a=mean(er);
end
%bounds for a suberset, contained in the identified set with probability
%1-alpha
G1hatbounds=boundsoriginal - epsilon*stdboundssigned


boundsoriginal


%%
%generating the graphs of figure 3 in the paper
figure('Visible','on',...
      'PaperPosition',[0 0 5 5],...
      'PaperSize',[5 5])

rectangle('Position',[-1,-1,2,2])
line([-1 1], [0,0],'Color',[0 0 0]);
line([0,0], [-1 1],'Color',[0 0 0]);

Ghat=rectangle('Position',[boundsoriginal(3),boundsoriginal(4),boundsoriginal(1)-boundsoriginal(3),boundsoriginal(2)-boundsoriginal(4)])
set(Ghat,'LineWidth',4,'EdgeColor',[0 0 0]);
text(boundsoriginal(3)+.03,boundsoriginal(4)+.07,'$\widehat{G}$','Interpreter','latex', 'Fontsize', 14)


G2hat=rectangle('Position',[G2hatbounds(3),G2hatbounds(4),G2hatbounds(1)-G2hatbounds(3),G2hatbounds(2)-G2hatbounds(4)])
set(G2hat,'LineWidth',3,'EdgeColor',[.8 .8 .8]);
text(G2hatbounds(3)+.03,G2hatbounds(4)+.07,'$\widehat{G}_2$','Interpreter','latex', 'Fontsize', 14)
%rectangle('Position',[G1hatbounds(3),G1hatbounds(4),G1hatbounds(1)-G1hatbounds(3),G1hatbounds(2)-G1hatbounds(4)])
 
 
slope=boundsoriginal(4)/boundsoriginal(1);
coneline;
set(cl,'LineWidth',4,'Color',[0 0 0]);
slope=boundsoriginal(2)/boundsoriginal(3);
coneline;
set(cl,'LineWidth',4,'Color',[0 0 0]);

slope=G2hatbounds(4)/G2hatbounds(1);
coneline;
set(cl,'LineWidth',3,'Color',[.8 .8 .8]);
slope=G2hatbounds(2)/G2hatbounds(3);
coneline;
set(cl,'LineWidth',3,'Color',[.8 .8 .8]);
 
set(gca, 'XTick', [-1:.2:1])
print(gcf, '-dpdf','Ghatmath05.pdf')
